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We study improved variants of Luscher's algorithm, including the cases of staggered fermions and finite density. 
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1. INTRODUCTION 

Simulation of full QCD, though essential, re- 
mains very expensive because of the large number 
of operations required by standard algorithms like 
Hybrid Monte Carlo (HMC) [1]. Cheaper algo- 
rithms are needed. Their cost should scale better 
than that of HMC as the volume of the system 
grows and the quark mass decreases. The algo- 
rithm proposed by Liischer [2] , which purports to 
be a controlled approximation to full QCD, repre- 
sents a viable alternative to HMC [3,4]. However 
it can be simply modified into more efficient, ex- 
act variants [5]. 

We summarize here our work on these variants, 
and extend the formulation to staggered fermions 
and finite density QCD. 

2. LUSCHER'S METHOD 

The method [2,3] approximates the fermionic 
determinant of 2 degenerate quark flavors as 



detQ 2 
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detP(Q 2 ) 



(1) 



where P(x) is a polynomial of even degree n ap- 
proximating 1/x in the interval (0,1]. Knowing 
the polynomial roots z k , k = 1, ..,n one can write 

n/2 

detP(Q 2 ) = const x JJ det(Q 2 -z k ) ] {Q 2 -z k )(2) 

k=l 

The inverse of each determinant on the r.h.s. 
can be expressed as a Gaussian integral over a 



bosonic field. We denote the Dirac matrix by D = 
(1 — kM), where M is the lattice Wilson hopping 
operator and k the hopping parameter. Liischer 
chooses for Q the hermitian operator Q = Cq Q, 
where Q = 75Z) and cq = (1 + 8/c) -1 . P(x) is 
constructed from Chebyshev polynomials and de- 
pends on an adjustable parameter e € (0, 1]. Its 
error is bounded uniformly in the segment [e, 1] 
[3]: 

H-^)!^^) (3) 

We assess the error of this method by measuring 
the fluctuations of det Q 2 P(Q 2 ) around 1, by cal- 
culating the eigenvalues of Q 2 over a sample of 4 4 
configurations [5]. 

3. EVEN-ODD SPLITTING 

By even-odd splitting the lattice sites, we fac- 
torize detQ 2 into equal, even and odd factors [5]: 



detQ 2 = det cg(l - K 2 M 2 ) 2 even 



(4) 
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with Co = (1+64k 2 ) -1 . It saves a factor 2 in mem- 
ory, but a bosonic update requires about as many 
operations as before. Nonetheless the smallest 
eigenvalue increases as seen from Fig.l. In Fig. 2 
we show the error decreasing faster with n after 
even-odd splitting. The gain is a factor 2 to 3. 

One can achieve further gains by rescaling Q 
to Q/cm- In [3] is assumed cm > 1 to guarantee 
that the spectrum is bounded by one. But for 
P finite one can take cm < 1, and the smallest 
eigenvalue of Q 2 will increase, leading to better 
convergence. Thus in Fig. 2 we tuned cm to 0.6 
instead of its default value of 1. 
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Figure 1. Left part: a typical spectrum of the 
Dirac matrix D in the complex plane (4 4 lattice, 
/3 = 6, k = 0.14). Right part: the spectrum is 
shown after even-odd preconditioning. 
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Figure 2. The magnitude of the error on the de- 
terminant as a function of n, in the original for- 
mulation (dotted line) and after even-odd precon- 
ditioning (solid line). 
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Figure 3. The magnitude of the error of the poly- 
nomial approximation P(x) as x varies from to 
1. The three approximations shown are the origi- 
nal hermitian approximation of Liischer, the non- 
hermitian approximation inside a circle as appro- 
priate for (3 = 0, and the non-hermitian approxi- 
mation inside an ellipse of aspect ratio 2, as ap- 
propriate for (3 = oo. All cases correspond to the 
same quark mass. 



4. NON-HERMITIAN VARIANT 

We look for an approximation detD « 
l/detP(D) where now P(z) of degree even is de- 
fined in the complex plane. Applying (2) for D 
and using the fact that D = 75.0^75 and 7! = 1 
one gets 

n/2 

detP(D) = const x JJ det(D - z k ) ] (D - z k ) (5) 

This approximation breaks down for detD nega- 
tive, which can only occur for very small quark 
masses. By this non-hermitian formulation one 
can simulate an odd number of quark flavors (or 
any number of non-degenerate flavors) . 

If the complex spectrum of D is bounded by 
the ellipse centered at (d, 0) , with large semiaxis 
a, and focal distance c, and if we construct P(z) 



by Chebyshev polynomials, we get [5] 

\ rt+l 

a + v a 2 — c 2 \ 

for all z within the ellipse defined above. We show 
in Fig. 3 the magnitude of the error along the real 
axis for n = 20 and the quark mass = 0.1. In 
the non-hermitian case the error falls off expo- 
nentially until it becomes uniformly oscillating in 
some interval. This is in sharp contrast with the 
hermitian approximation of Liischer. 

We measure the fluctuations of (DP(D)) 2 by 
calculating the eigenvalues of D applying the 
same methodology as before. In Fig. 4 we show 
that the gain over the hermitian case is up to a 
factor 4 to 5 at (3 = 0. Note that simulating one 
flavor to the same accuracy would require half as 
many bosonic fields n. 
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Figure 4. The magnitude of the error on the de- 
terminant as a function of n, in the original for- 
mulation (dotted line) and in the non-hermitian 
variant (solid line). 



5. METROPOLIS TEST 

Liischer's original proposal includes the moni- 
toring of the error, and the possibility of obtain- 
ing exact results by re-weighting the Monte Carlo 
measurements of each observable. It is more at- 
tractive however to correct for the error through 
a Metropolis test, just like HMC corrects for the 
error of the guiding Hamiltonian. 

A first attempt along these lines [6,4] used the 
Lanczos algorithm to obtain eigenvalues of Q 2 . 
The error det Q 2 P(Q 2 ) is then easy to calcu- 
late. However the cost of this procedure grows 
prohibitively like the square of the volume V of 
the lattice, unless one accepts an "error on the 
error" . 

We proposed instead [5] to use a noisy, unbi- 
ased estimator of the ratio of errors between the 
new and the old configurations. The ratio to esti- 
mate is det(D' P(D')) 2 /det(DP(D)) 2 , calling D' 
and D the Dirac operators for the new and old 
configurations respectively. Taking the denomi- 
nator as a partition function, this ratio can be 
rewritten 



< e -r,\w^w-i)r, > 



(7) 



where the average <> is taken over all Gaussian 
vectors jj, and W = [D 1 P (D')] 1 D P (D) . It is 
sufficient to estimate this ratio by taking one 



Gaussian rj only at the end of each "trajectory" , 
thus solving one linear system as in HMC. It is 
easy to show that the associated cost, expressed 
in number of update sweeps, is constant as the 
volume of the system grows or the quark mass 
decreases. The proof of detailed balance is more 
subtle. 

Consider a generic update U — > U' of the gauge 
fields, consisting of a "trajectory" of Liischer up- 
dates followed by a Metropolis step. Detailed bal- 
ance is satisfied if 

Pl(U -> U 1 ) Pacc(U'\U) = e- s °( u 'Uet 2 (D>) 
P L (U'^U) P acc (U\U') e- s °( u )det 2 (D) U 

where Pl and P acc are the probabilities for the 
set of Liischer steps and the Metropolis step. Let 
us assume that the Liischer updates satisfy de- 
tailed balance with respect to the Liischer action 
e~ 9 det~ 2 (P '(D)) . This can be simply achieved 
by taking care that the order of the successive U 
and 4>k updates in a trajectory can be reversed. 
Then eq.(8) implies 



Pacc(U'\U) det 2 (D'P(D')) 



(9) 



Pacc(U\U>) det 2 (DP(D)) 
If one choses one Gaussian vector r\ and simply 
takes P acc (U'\U) = mm(l,e-" t ( w ' tw '- 1 )'>), then 
P acc (U\U') = mmtLe-"^^ 1 ^ 1 )- 1 )"), and 
detailed balance will not be satisfied. 

However it is simple to restore detailed balance 
with the following prescription: 
. if, say, S g ([/')> S g (U), 
P acc (U'\U) =min(l,e- r > i (w^w-i)^ 
• otherwise, 

Pacc(U'\U) = 



l ( 1; e+» t ((w- 1 tw- 1 )-i)')) 



6. STAGGERED FERMIONS 

Liischer's method can be applied to staggered 
fermions. In this case the Dirac matrix has the 
form D = ml + iB, where m is the bare mass and 
B is a hermitian matrix with extreme eigenvalues 
±A. Then, in the hermitian formulation we have 
Q 2 = D^D = m 2 l + B 2 . Using (6) with a = c = 
X 2 /2 and d = m 2 + X 2 /2 the error bound will be 

, , \2(n+l) 

\l-zP(z)\<2 - —= ) (10) 



2 + y/m 2 + X 2 
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for z e [m' 2 ,m 2 + A 2 ]. 

Similary as for Wilson fermions, one can also 
define a non-hermitian approximation using (2) 
for D and the identities D = E 2 = 1, 

where £ is ±1 on even/odd sites. In this case, 
taking a = 0, c = i\, and d = m we get the 
same error bound as above, but now for z G [in — 
i\, m + 

Thus, there are no savings from a non- 
hermitian approximation. The situation is the 
same as for the calculation of staggered fermion 
propagators, where no savings over CG can be 
achieved by a non-hermitian solver [7]. 

7. FINITE DENSITY QCD 

There are several difficulties to simulate finite 
density QCD: i) since the determinant is not real 
anymore, its representation as a Gaussian inte- 
gral over bosonic fields may not be convergent. 
Additional problems are: ii) how to get the phase 
of the determinant without calculating it directly 
and iii) how to deal with the large MC fluctu- 
ations caused by the complex phase (sign prob- 
lem). We address here only problems i) and ii). 

If we express exp(n) = cosh(n) + sinh(/i), then 
we can decompose the quark matrix as D = 
Di + D 2 , where D 2 contains the sinh(/i) term 
and D\ = 75D175, D\ = — 75-D275. If we assume 
an even degree polynomial P(z) defined in the 
complex plane and apply (2) for Di + D2, we get 

detP(D) = 

n/2 

const x JJ det(D 1 - D 2 - z k )^(D 1 + D 2 - z k ) = 

k=l 

n/2 

const x JJ dei[(L»i - z k )\D 1 - z k ) - D\D2 + 

k=\ 

(D 1 -z k )^D 2 -Dl(D 1 -z k )} (11) 

In the last expression (Di —z k )^(D\ — z k ) — D\Di 
is hermitian, whereas {D\ — z k )^ D 2 — d\{D\ — z k ) 
is antihermitian. There exists an interval of n 
around for which the hermitian part is positive 



and thus the corresponding bosonic Gaussian in- 
tegral is convergent. In this case we can esti- 
mate the magnitude of the determinant from the 
hermitian part of (11), whereas the antihermitian 
part will give an estimation of its phase. This way 
we can approximate the determinant without cal- 
culating it directly. 
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